setwd("C:/Users/dzhao1/Desktop/LMM")
cnames=c("id","trt","y0","y2","y4","y6","y8","y10","y12")
exercise.wide=read.table("stren.txt", col.names=cnames, na.strings = ".")
#reshape the data from wide format to long format
exercise.long <- reshape(exercise.wide, idvar="id", varying=c("y0","y2","y4","y6","y8","y10","y12"), 
                        v.names="y", timevar="time",time=c(0,2,4,6,8,10,12), direction="long")

attach(exercise.long)
mean.na.rm=function(x) mean(x,na.rm=T)
# Spaghetti Plot of individiual trajectories
interaction.plot(time, id, y, type="b", fun=mean.na.rm, 
                 xlab="Time (in days)", ylab="Boday Strength", 
                 main="Spaghetti Plot", col=c(1:37), legend=F)
#mean profile plot
interaction.plot(time, trt, y, type="b", fun=mean.na.rm,
                 xlab="Time (in days)", ylab="Boday Strength", 
                 main="Mean Response Plot")
#######################################################
#install package
install.packages("lme4") 
#load package
library("lme4")
#Fit a LMM with fixed time, trt, and random patient,REML estimation
excercise1 <- lmer(y ~ time + trt + (time | id), exercise.long, REML=T)
#Summary of the results (fixed effects, AIC, loglik, etc)
summary(excercise1)
#Get 95% CI for fixed effects and variances
confint(excercise1)
#estimated random effects
ranef(excercise1)
#individual effect
coef(excercise1)
#individual prediction (fixed effects + random effects)
predict(excercise1)
#Mean prediction (fixed effects, without random effects)
predict(excercise1, re.form=NA)
#Residual plot
plot(excercise1)
qqnorm(resid(excercise1))
qqline(resid(excercise1))